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On  the  Iterative  Boundary  Element  Method 


B  D.  Cahan  and  O.  Lafe 

Case  Center  for  Electrochemical  Sciences,  Case  Western  Reserve  University, 
Cleveland,  Ohio  44106. 

ABSTRACT 

A  generalized  iterative  scheme  is  developed  for  the  solution  of  a  large  class  of 
boundary  integral  equations.  The  code  is  an  excellent  method  for  solving  prob¬ 
lems  with  huge  coefficient  matrices.  The  scheme  avoids  the  use  of  large  computer 
storage  spaces  and  a  direct  inversion  of  the  matrices  is  also  not  required.  The 
result  is  a  considerable  improvement  over  existing  boundary  element  schemes 
in  terms  of  computer  running  time  and  accuracy  especially  for  problems  with 
complex  shapes  and  non-linear  boundary  conditions. 

INTRODUCTION 

The  boundary  element  method  (BEM)  has  evolved  within  the  last  two  decades 
as  a  powerful  tool  for  solving  a  variety  of  problems  in  physics  and  engineering. 
The  automatic  reduction  in  the  dimensionality  of  the  problem  is  one  of  the  major 
advantages  BEM  has  enjoyed  over  the  domain  techniques  such  as  finite  element 
and  finite  difference  methods 

One  key  difficulty  in  the  use  of  BEM  for  serious  real-life  cases  is  the  problem  of 
ill-conditioning  in  the  global  matrix  when  a  fine  resolution  is  required  in  certain 
localities  of  the  domain  boundary.  Examples  include  the  formation  of  fractals 
in  electrodeposition  problems,  cracks  in  fracture  mechanics  and  fingering  in  soil 
moisture  studies.  Zoning  has  been  suggested  as  a  way  of  avoiding  such  problems 
(see  Lafe  et  a l  [3|).  However  zoning  increases  the  complexity  and  computational 
effort  of  the  boundary  solution  process  since  artificial  boundaries  will  then  have 
to  be  introduced  in  the  problem  domain.  The  effect  is  an  impairment  of  the  user- 
friendliness  in  terms  of  the  ease  of  data  preparation  normally  associated  with 
BEM  codes. 

Furthermore  the  incorporation  of  non-linear  mixed  boundary  conditions  (e  g. 
kinematic  conditions  in  free-surface  fluids  problems  or  Tafel  relationships  in  elec¬ 
trochemistry)  is  a  rather  difficult  task  in  BEM  The  solution  process  in  such 
cases  most  often  requires  the  inversion  of  the  coefficient  matrix  several  times. 
Such  clumsy  schemes  destroy  the  traditional  advantages  of  BEM. 

In  this  paper  we  propose  a  simple  but  elegant  solution  scheme  suitable  for  a 
variety  of  boundary  integral  problems.  We  call  it  the  Iterative  Boundary  Element 
Method  (I-BEM)  in  that  the  entire  solution  process  is  iterative  even  for  linear 
equations  with  linear  boundary  conditions.  The  kernel  of  the  method  was  first 
suggested  by  Cahan  jll  and  then  implemented  by  Cahan  et  al  \2\.  The  origi¬ 
nal  scheme  involved  locating  integration  points  near  the  boundary  just  inside  of 
the  problem  domain  and  the  use  of  finite  difference  as  a  predictor  in  the  itera¬ 
tive  prc  ess.  The  scheme  is  generalized  in  this  paper  to  accommodate  arbitrary 
choices  of  integration  points  on  or  off  the  boundary. 

/)-/• 


A  comprehensive  description  of  the  method  is  given  here  and  we  show  how  a 
relaxation  factor  A  can  be  effectively  used  to  accelerate  the  rate  of  convergence. 
We  also  examine  the  relative  performances  of  I-BEM  and  BEM  by  solving  a 
problem  connected  with  the  electrodeposition  at  an  external  corner.  This  problem 
is  characterized  by  exponentially  strong  boundary  conditions. 

GOVERNING  DIFFERENTIAL  EQUATION 

For  a  domain  H  consisting  of  the  boundary  T  we  look  at  a  general  class  of  problems 
governed  by  the  differential  equation: 

£W(x.t)|  =  -r,(x'.  t')  S(xr,  x")  (1) 

where  £'  is  some  linear  partial  differential  operator  in  both  space  (x)  and  time  (t), 

<t>'  is  the  dependent  variable  (e  g.  voltage,  pressure,  temperature,  displacement, 
potential,  stream  function  etc),  6  is  the  Dirac  delta  function  while  7'  represents 
the  strengths  of  point  actions  or  sources  in  Q. 

By  taking  initial  conditions  as  reference  values  we  define  4>  —  fn°  4>'  exp (-st)dt 
where  s  is  the  Laplace  transform  parameter  and  use  these  to  transform  equation 
(l)  into: 

£-o(x.  s)j  =  7(x.  s)  <5(x',  x")  (2) 

in  which  £  is  now  a  differential  operator  in  space  while  7  is  the  transform  of  7' . 

In  subsequent  statements  we  drop  the  suffix  s  without  any  loss  of  generality. 

The  conditions  associated  with  the  boundary  F  can  in  general  be  written  in 
the  form: 

f(O.Cld>  i.x)-0  (3) 

where  £  is  another  differential  operator  and  /  is  some  prescribed  function.  On 
a  Di  rich  let  boundary  /  =  f{o,x)  while  on  a  Neuman  segment  /  =  f(P\o  j,x)  at 
most. 

INTEGRAL  EQUATION 

The  first  task  in  a  BEM  process  is  to  convert  the  governing  partial  differential 
equation(l)  into  a  suitable  integral  equation.  This  is  achieved  by  multiplying 
the  equation  by  a  function  G,  integrating  the  ensuing  expression  over  0  and 
invoking  the  Green’s  identities  or  the  Betti's  reciprocal  theorem.  The  result  of 
such  manipulations  is  the  integral  equation: 

j  d>(x)£’(G(x,x')|  dU  ~  J  (d>(x')G’(x,  x') -<2>'(x')G(x,  x')  dr G(x.x")7(x") 

(4) 

where  £*  is  the  adjoint  to  operator  £  while  G‘  and  o'  are  derivative  functions 
of  G  and  4>  respectively.  If  we  choose  the  function  G  such  that: 

£'  G(x.  x')j  —  A(x.  x')  (5) 

then  equation(d)  degenerates  into 

od>(x)  =  j  (4>{x')G' (x.  x')  -  o'(x')G(x.  x')  dr ^  G(x,  x").r(x")  (6) 
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where  a  is  the  Cauchy  principal  value  of  the  integration  of  the  Green’s  function 
singularity. 

COEFFICIENT  MATRICES 

To  use  the  boundary  element  method  the  boundary  F  is  subdivided  into  a  finite 
number  of  elements  and  suitable  interpolation  functions  are  chosen  to  represent 
the  distribution  of  the  dependent  variables  on  the  boundary.  For  function  <t>  on 
T  we  write: 

d>(x)  =  H  Nt(*)4>i  (7) 

I 

where  <£>,  are  the  values  of  <t>  at  the  discrete  points  on  the  boundary  and  Nx 
are  shape  functions.  A  similar  expression  can  be  written  for  <t>' .  To  perform 
boundary  integrations  we  select  a  finite  number  of  points  to  serve  as  integration 
origins  such  as  x  in  equation(6).  In  most  BEM  implementations  these  points 
are  located  along  the  boundary  segments  and  usually  fall  on  the  nodes  formed 
at  element  intersections.  In  a  few  implementations  (the  so-called  non-singular 
formulation)  the  integration  points  are  selected  outside  H.  In  1-BEM  x  can  be 
placed  on  T.  inside  or  outside  of  H.  The  number  of  independent  integration  points 
selected  is  the  same  as  the  number,  M.  of  all  unknown  4>  or  <p‘  on  the  boundary. 
When  these  integrations  are  performed  equation(4)  becomes: 

Vf  M 

]La‘A  =  +  c*  ‘  =  1. 2,  •  •  • ,  Af 

;  - 1  ;  =  i 

where 

atJ  =  /  ,Vt'(x')C(x.x')  dr 

b,}  =  f  N,(x')G’(x.x')  dr 

c,  =  Yi  -r(xfc)G(x,xfc) 

k 

In  a  conventional  BEM  the  boundary  conditions  will  be  introduced  into  the 
system  of  equations(8)  and  then  assembled  into  the  form: 

M 

Y,  d*;u;  ~  «  =  12,  •••.A/  (9) 

where  u,  represents  the  unknown  <t>  or  <t>'  at  the  M  nodes.  Equation(9)  is 
then  inverted  to  obtain  these  unknown  quantities.  For  problems  with  non-linear 
boundary  condition?  equation!!})  becomes  essentially  a  non-linear  system  of  equa¬ 
tions  and  the  assembling  process  may  even  have  to  be  repeated  several  times  in 
order  to  take  full  advantage  of  available  linear  matrix  inversion  routines.  The 
approach  in  l-BEM  is  different  since  the  assembling  symbolized  by  equation(9) 
is  never  carried  out  Rather  equalion(8)  is  utilized  directly. 


(8) 
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I-BEM 

To  solve  for  the  unknown  <t>  and  <f>'  at  the  nodes  we  start  by  guessing  the  values. 
4>t  and  for  these  quantities.  (Zeroes  may  be  used  for  unfamiliar  problems.  The 
additional  number  of  iterations  required  is  usually  small  (2-3)  when  compared  to 
the  case  when  a  more  ‘intelligent’  initial  guess  is  made.)  For  Dirichlet  (exact  <i>,) 
and  Neuman  (exact  4>\ )  type  problems  the  known  quantities  can  be  used  onre 
and  for  all  in  equation(8)  and  the  results  absorbed  into  the  coefficients  c, .  For 
linear  mixed  boundary  conditions  one  will  normally  eliminate  either  <p  or  <£>'  for 
the  other.  In  non-linear  mixed  conditions  the  iterative  process  readily  accepts 
the  incorporation  of  any  root-finding  routine  in  the  solution  process.  When  the 
guesses  are  used  in  the  equation  written  for  the  «-th  node  as  origin  of  integration 
[see  equation(8)|  the  result  is  an  error: 

AY  AY 

s,  =  *■  c«  (10) 

;=•  ; -i 


Assuming  the  largest  contributor  t,o  e,  comes  from  the  j  =  i  terms  the  j  ^ 
t  terms  in  equation(IG)  can  be  suppressed  to  obtain  a  suitable  expression  for 
updating  the  unknown(s)  at  the  i~th  node.  The  ensuing  update  equations  are: 


Neuman  Boundary 

-  t 

<t>,  =  <t>.  -  A 

Ox, 

(11) 

Dirichlet  Boundary 

<t>\  -  O’  *  A-'1 

hxx 

(12) 

Mixed  Boundary 

o„(o,  -  -  Oft)  =  -As,  (13) 

where  A  is  an  over-relaxation  factor  whose  optimal  choire  is  examined  below.  For 
a  segment  with  a  mixed  boundary  condition  equation(  13)  is  to  be  solved  with  the 
relevant  prescribed  condition  as  shown  by  equation(3) 

The  above  process  is  repeated  for  all  points  until  a  suitable  convergence  cri¬ 
terion  is  satisfied. 


OPTIMAL  A- VALUES 

An  exact  analysis  of  the  optimal  A  is  difficult  to  carry  out  We  present  an  approx¬ 
imate  analysis  which  provides  a  sense  of  the  best  estimate  for  the  over-relaxation 
factor. 

We  write  error  equation!  10)  for  two  iteration  levels  k  and  Jtr 1  By  subtracting 
the  ensuing  expressions  the  result  is: 


.fc  +  i 


Af 


;  i 


o 


It- 1 

1 


O, 


AY 

;  =  t 


(14) 
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Using  equations  (11)  and  (12)  in  (14)  and  rearranging  we  obtain: 

fc+l  M 


For  convergence  we  must  enforce: 


J  =  ! 

'  I 


/f*\ 

_aj7  +  h} 

w; 

.ail  b}} 

(15) 


<  1 


Assuming  a  uniformity  in  segment  errors  (c*  *s  etfc)  equation  (15)  when  combined 
with  the  above  condition  produces: 


.vr 

j  =  i 


b 


n 


(16) 


Since  condition  (16)  would  provide  different  limitations  on  the  value  of  the  As 
for  different  segment  and  our  analysis  assumed  an  a  priori  constant  A  for  all 
segments  the  optimal  choice  has  to  be  the  value  that  satisfies  the  condition  for 
all  segments.  That  is  in  addition  to  being  positive  the  optimal  value  has  to  be 
such  that: 

l  /  r 


A  < 


' u 


t  = 


(17) 


Af 

■n  2 :  E 

Numerical  experiments  have  confirmed  the  validity  of  the  above  convergence  cri¬ 
terion.  The  result  shows  the  optimal  A  is  dependent  on  geometry  and  the  number 
of  boundary  elements 

NUMERICAL  RESULTS 

The  following  numerical  tests  have  been  performed  for  a  steady  state  boundary 
value  problem  governed  by  Laplace's  equation  («.c.  £  —  V2  and  ?  —  0).  The 
domain  of  study  is  as  depicted  in  Fig(l). 

The  problem  is  connected  with  the  electrodeposition  process  in  which  the 
condition  on  a  portion  the  boundary  is  governed  by  the  Butler- Volmer  equation: 

d<t> 
dn 

where  k,  tq  and  iq  are  coefficients  dependent  on  the  material  and  kinetics  of 
the  reaction.  In  the  special  case  when  \d>\  >  >  I  iq  the  Butler- Volmer  equation 
degenerates  to  the  so-called  Tafel  equation  In  all  the  simulations  linear  shape 
functions  have  been  used.  Error  plots  for  varying  values  of  relaxation  factors  and 
number  of  of  boundary  elements  are  shown  in  Fig(2) 

The  iterative  process  was  carried  out  over  3  passes.  The  maximum  observed 
error  is  4.9  x  10~V  The  minimum  is  2  I  ■  10~5  for  the  3  iterations  The  optimal 
relaxation  factor  generally  lies  between  0  6  -  0  75  The  higher  value  seems  to 
be  more  appropriate  for  problems  with  relatively  large  number  of  elements.  The 
convergence  rate  is  reduced  significantly  when  a  small  factor  <  0.2  is  used  for 
many  test  runs  that  were  made  The  0  75  upper  limit  is  rigid  for  most  problems 
as  severe  instabilities  are  common  for  relaxation  factors  just  slightly  higher  than 
this  value 


n\e 


,  -  I'J*! 
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Figure  2:  Error  Plots  for  various  A  values  and  element  sue.  Normalized  error 

=  (<_<min)/(fma:'(mm) - - - —  10  elements  segment; - * - -  80 

elements/ segment 


I-BEM  VERSUS  BEM 


In  using  BEM  the  presence  of  the  Tafel  segment  means  the  the  final  global  matrix 
is  nonlinear.  The  solution  process  then  even  in  that  case  still  has  to  be  iterative  as 
mentioned  earlier.  The  comparative  central  processing  unit  (cpu)  time  plots  for 
both  methods  are  shown  in  Fig(3).  Both  processes  were  limited  to  3  iterations. 
While  the  typical  error  in  I-BEM  (with  a  A  value  of  0.7)  was  of  the  order  JO'4 
that  of  BEM  was  of  order  10'2  at  the  end  of  the  3  iterations.  The  significant 
advantage  enjoyed  by  I-BEM  in  terms  of  computational  speed  is  quite  evident 
from  the  two  curves  in  Fig(3)  especially  as  the  number  of  elements  is  increased. 

CONCLUSIONS 

The  iterative  boundary  element  method  has  been  presented  in  this  paper  as  a 
powerful  tool  for  solving  a  wide  class  of  boundary  value  problems.  The  advan¬ 
tages  of  the  new  scheme  are  much  evident  in  large  systems.  A  given  problem 
with  M  unknowns  on  the  boundary  typically  requires  C>(A/3)  arithmetic  compu¬ 
tations  for  the  direct  matrix  inversion  used  in  conventional  BEM.  In  I-BEM  the 
computational  effort  is  C([m,lzXf2)  where  is  the  total  number  of  iterations 
required  to  reach  convergence.  Thus  I-BEM  enjoys  a  computational  advantage 
of  the  order  \f  fmai  over  PEM  In  general  the  number  of  iterations  necessary  is 
orders  of  magnitude  smaller  than  the  number  of  unknowns.  Lafe  and  Cheng  |4| 
have  reported  the  simulation  of  a  stochastic  flow  problem  in  which  the  computa¬ 
tional  time  using  l-BEM  for  1600  unknowns  was  5%  of  the  corresponding  time 
required  by  BEM 
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Figure  3:  Comparative  run-time  plots  for  I-BEM  and  BEM.  -  — 

BEM; - - - —  I-BEM 
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